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The  objective  of  thi«  work  has  been  to  develop  a  direct  numerical  simulation  capability  for  compressible, 
non-eanonical  wall-bounded  flows  which  is  directly  applicable  to  actual  flight  vehicles.  Data  from  the  simulations 
is  used  to  educe  information  about  the  organised  motions  in  the  flow  and  how  they  axe  affected  by  the  extra 
strain  rates  due  to  concave  wall  curvature.  A  second-order  finite  volume  approach  and  both  fourth-order  and 
sixth-order  compact  differences  are  available  for  the  spatial  derivatives.  The  algorithm  allows  for  generalised 
coordinates  so  that  simulations  about  complex  geometries  can  be  performed.  Extensive  testing  was  done  for  two 
subgrid-scale  models,  the  compressible  Smagorinsky  and  strneture  function  models,  for  a  supersonic  boundary 
layer.  As  the  subgrid-scale  viscosity  also  behaves  as  an  artificial  viscosity  for  central  differencing  codes,  the 
study  reinforced  the  need  for  high-order  dissipation  models  in  the  simulation  code.  Validation  of  the  numerical 
method  was  performed  for  flow  over  a  concave  surface  for  subsonic  flow  in  which  streamwise  vortices  develop  in 
the  boundary-layer  due  to  the  Gortler  instability.  The  geometry  and  flow  conditions  closely  approximated  the 
incompressible  experiments  of  Swearingen  and  Blackwelder.  The  simulations  captured  the  essential  features  of 
the  experiments  in  which  counter-rotating  vortices  developed  near  the  wall.  As  these  streamwise  vortices  develop, 
the  laminar  boundary  layer  on  the  concave  wall  rapidly  becomes  three-dimensional.  Higher-momentum  fluid  is 
pulled  towards  the  wall  from  the  outer  flow  as  vortices  grow  downstream.  The  low-speed  fluid  lying  between  the 
vortices  induces  low-momentum  fluid  away  from  the  wall,  leading  to  inflectional  streamwise  velocity  profiles,  in 
good  agreement  with  the  experiments.  X)TIC  QUALITY  INSPECTED  3 
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1  OBJECTIVES 

The  objective  of  this  research  effort  was  to  develop  direct  numerical  simulation  (DNS) 
and  large-eddy  simulation  (LES)  techniques  to  simulate  compressible,  non-canonical  wall- 
bounded  flows  which  are  more  directly  applicable  to  actual  flight  vehicles.  The  data  from  the 
simulations  is  used  to  educe  information  about  the  organized  motions  in  the  flow  and  how 
they  are  affected  by  the  extra  strain  rates  due  to  the  wall  curvature.  Also,  the  data  from  the 
simulations  can  be  used  for  improvements  to  the  Reynolds  stress  transport  equations  and  the 
two-equation  k  —  e  turbulence  models  for  compressibility  and  curvature  effects.  Specifically, 
the  research  objectives  are  the  following: 

•  Develop  a  high-order  compact  difference  algorithm  for  solution  of  the  compressible, 
Navier- Stokes  and  energy  equations  for  simulation  of  spatially  evolving  three-dimensional 
turbulent  flow. 

•  Investigate  boundary  condition  compatibility  issues  and  dissipation  models  for  the 
higher-order  differencing. 

•  Validate  the  code  for  small  amplitude  disturbances. 

•  Evaluate  subgrid-scale  models  for  wall-bounded  compressible  flows. 

•  Compute  the  turbulent  flow  over  a  flat  plate  at  a  supersonic  Mach  number  and  compare 
with  incompressible  simulations. 

•  Compute  the  transitional  flow  over  a  concave  surface. 

•  Characterize  the  nonlinear  stages  of  transition  to  determine  the  effects  of  the  extra 
strain  rates  applied  by  the  curved  surface. 
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2  GOVERNING  EQUATIONS 
2.1  Governing  Equations 

The  compressible  transitional  and  turbulent  flow  of  an  ideal  gas  is  governed  by  the  conti¬ 
nuity,  momentum  and  energy  equations.  The  strongly  conservative  form  of  these  governing 
equations  will  be  solved.  The  equations  can  be  written  as  follows,  where  uk  =  [u,  v,  w]T 
is  the  velocity  and  p,  p,  T,  /x,  k,  and  Cv  are  the  density,  pressure,  temperature,  viscosity, 
thermal  conductivity,  and  specific  heat  at  constant  volume  respectively: 


dp  (dpuk) 

at  dxk 

(1) 

d(puk)  d(pukui)  _  dp  derki 

dt  dxi  dxk  dxi 

(2) 

8(PT)  8(p»tT)  8u„  8  (  BT\ 

Cv  dt  +Cv  dxk  -  Pdxk  +  dxk  [Kdx  J  +  * 

(3) 

The  viscous  stress  <7 ki  and  the  viscous  dissipation  $  are  defined  by 

2  duj  ( duk  dui  \ 

at'~~Tar,s‘,+'i{a^,  +  ari) 

(4) 

2  ( duj\ 2  ( duk  dui\  duk 

r[aj  +',U.+&J&r 

(5) 

The  equation  of  state  is 

p  =  pRT. 

(6) 

The  viscosity  and  thermal  conductivity  are  varied  with  temperature  using  Sutherland’s  law 
for  the  viscosity  and  a  fourth-order  polynomial  curve  fit  to  experimental  data  for  the  thermal 
conductivity.  The  specific  heat  at  constant  volume  is  assumed  constant.  The  physical  domain 
(x,  y,  z)  is  mapped  to  the  computational  domain  (£,  77,  £)  to  allow  for  arbitrarily  shaped 
bodies. 

2.2  Favre  Filtered  Equations 

In  LES  the  flow  variables  are  decomposed  into  a  large-scale  component  and  a  subgrid-scale 
component.  The  large-scale  component  is  defined  by  a  filtering  operation: 
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where  Gi  is  a  filter  function  in  the  *th  direction  and  D  is  the  domain  of  the  fluid.  Two 
commonly  used  filters  are  the  sharp  Fourier  cutoff  filter,  defined  in  physical  space  as 

G,(x,-,x-)  =  2  sin[jr(x,  -  x{)/At] /x(x,- -  x-)  (»  =  1,3)  (8) 

and  the  Gaussian  filter  defined  as 

Gi(xi,  x'i)  =  (6/* A,)1/2  exp  [-6  (x,-  -  x')2  /A?]  (*  =  1,3),  (9) 

where  A ;  is  proportional  to  the  grid  size  in  the  tth  direction.  For  the  Fourier  cutoff  filter 

the  subgrid  scale  field  contains  the  velocity  due  to  all  the  structures  with  wavenumber  |Jfe,|  > 
7t/A,.  The  Gaussian  filter  is  a  more  global  filter  since  a  wider  range  of  scales  contribute  to 
the  subgrid-scale  velocity. 

Favre  filtering  is  used  to  decompose  the  turbulent  field:  T  =  T -\-T' ,  where  T'  is  the  SGS 
part  of  T  and  the  Favre  filter  is  defined  by  T  —  pT /p.  The  filtering  operation  is  applied  to 
the  governing  equations  which  yields  filtered  equations  for  the  large  eddies: 


dp  .  d(puk) 
dt  dxk 

(10) 

d(puk)  d(pukui)  dp  d<Tki  dm 

dt  dxi  dxt  dxi  dxi 

(11) 

r9(pT) 

Cv  dt 

Sifmd)  Bu„  -  B  (.BT\  r8Qk 

+  C-  dxk  C°8xk 

(12) 

p  =  pRT . 

(13) 

The  SGS  stress  tensor  m  and  the  SGS  heat  flux  Qk  are  defined  by 

Tkl  =  —p{ukUl  —  XLkU l  +  tijfcflf  -f  u\uk  +  u'kU'l) 

(14) 

Qk  =  p{ukT  -u kT  +  j~f  +  ^jf'  +  ^T'). 

(15) 

2.3  Subgrid— Scale  Modeling 

The  capabilities  of  two  SGS  stress  models  to  accurately  model  wall-bounded,  compressible, 
turbulent  flows  are  investigated.  The  SEZHu  model  and  the  structure  function  model  are 
chosen  due  to  their  success  in  modeling  the  SGS  physics  of  compressible  isotropic  turbulence. 
A  brief  description  of  these  models  follows. 

The  SEZHu  SGS  model  was  derived  by  Speziale  et  a 1,  (1988)  for  compressible  isotropic 
turbulence.  This  model  is  based  on  the  Favre-filtered  equations  of  motion  of  an  ideal  gas. 
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Piomelli  et  al.  (1988)  has  shown  that  the  structure  of  the  subgrid  scales  depends  strongly  on 
the  type  of  filter  used  and  consistency  between  the  model  and  filter  is  necessary  for  accurate 
results.  Following  the  work  of  Piomelli  et  al.  (1988),  only  the  Smagorinsky  portion  of  the 

SEZHu  model  is  used  with  the  Fourier  cutoff  filter.  Hence,  the  SGS  stress  model  is  of  the 

form 

rw  =  (l  -  e-y+'25)2  2CRpA2Ill/2  (&,  -  (16) 

and  the  SGS  heat  flux  is  given  by 

=  (17) 

where  Ski  =  (duk/dxi  -f  dui/dxk)  /2  is  the  Favre  filtered  rate  of  strain  tensor  and  11$  = 
Smn Smn  is  its  second  invariant.  Cr  is  the  compressible  Smagorinsky  constant,  Pr t  is  the 
turbulent  Prandtl  number,  and  A  =  (AxAyAz)1/3.  The  first  term  in  these  equations  rep¬ 
resents  a  Van  Driest  wall  damping  and  y+  is  the  distance  from  the  wall  nondimensionalized 
by  wall  shear  velocity  and  kinematic  viscosity. 

The  second  SGS  model  is  the  structure  function  model  (Comte  et  al.  1990)  which  is  based 
on  the  concept  of  spectral  eddy  viscosity  and  spectral  eddy  diffusivity.  A  local  kinetic  energy 
spectrum  is  calculated  in  terms  of  the  local  second-order  velocity  structure  function.  The 
structure  function  model  is  of  the  following  form  for  the  SGS  shear  stress  and  heat  flux: 

ru  =  (l  -  e-"+/25)2  CRpAyfF2(x,Ax,t)  x  (§u  -  ^SmmSk^  (18) 


where 

7 7  =  (I  W*,<)  -  <Ki+ '.0IP)„r1|=iI  •  (20) 

Again,  Van  Driest  wall  damping  is  applied. 

3  NUMERICAL  APPROACH 

A  large  extent  of  this  research  effort  has  been  directed  towards  development  of  the  de¬ 
velopment  of  the  high-order  compact  difference  simulation  code  for  generalized  curvilinear 
coordinates,  as  outlined  in  the  objectives.  Towards  this  end,  several  studies  were  performed 
to  determine  the  most  appropriate  numerical  techniques.  The  numerical  algorithm  and  sev¬ 
eral  key  numerical  issues  are  highlighted  below,  including  application  to  specific  test  cases. 
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For  the  time  integration  scheme,  the  second-order  Runge-Kutta  methods  by  Jameson 
(1980),  the  low  storage  third-order  Runge-Kutta  methods  by  Wray  (1991)  and  Williamson 
(1980),  and  the  classical  fourth-order  Runge-Kutta  method  have  been  evaluated.  The  third- 
and  fourth-order  schemes  were  tested  for  a  simple  ordinary  differential  equation  and  the  pre¬ 
dicted  third-  or  fourth-order  accuracy  was  achieved.  The  third-order  Runge-Kutta  method 
requires  fewer  computations  and  less  storage  than  the  fourth-order  Runge-Kutta  method 
and  is  being  used  in  the  full  simulations. 

For  the  spatial  differencing  procedure,  a  second-order  finite  volume  approach  and  fourth- 
order  or  sixth-order  compact  differences  are  available.  A  generalized  coordinate  system  is 
being  used  to  allow  for  arbitrarily  shaped  bodies.  A  cell-centered  scheme  is  employed  for 
the  second-order  finite  volume  method.  The  compact  difference  scheme  is  applied  at  the 
nodes,  not  the  cell  centers,  since  higher  order  interpolation  schemes  would  be  required  for  a 
cell-centered  scheme,  resulting  in  many  additional  calculations. 

Both  the  fourth-order  and  sixth-order  compact  differences  form  a  tridiagonal  system  of 
equations.  The  difference  approximations  for  the  near  boundary  nodes  are  necessarily  one¬ 
sided  and  consistent  with  the  overall  compact  difference  scheme.  Following  the  work  of 
Lele  (1990),  Carpenter  (1991)  and  Pruett(1993),  the  following  schemes  are  employed  at  the 
boundaries.  For  the  fourth-order  compact  differences,  a  3-4-3  stencil  is  used,  which  repre¬ 
sents  fourth-order  differencing  in  the  interior  of  the  mesh  and  third-order  at  the  boundaries. 
For  the  sixth-order  compact  differences,  either  a  3, 4-6-4, 3  stencil  or  a  5, 5-6-5, 5  stencil  is 
used  in  which  sixth  order  compact  differences  are  used  in  the  interior.  For  the  first  method, 
third-order  compact  differences  are  used  at  the  boundaries  and  fourth-order  differences  are 
used  at  the  first  interior  points.  Similarly,  the  second  sixth-order  method  employs  fifth- 
order  differences  at  the  boundaries  and  the  first  interior  points.  The  code  is  programmed  so 
that  either  set  of  boundary  conditions  can  be  used  for  the  sixth-order  compact  differences. 
However,  Pruett  has  shown  that  stability  is  maintained  when  the  5, 5-6-5, 5  scheme  is  used 
in  the  streamwise  direction  and  the  3,4— 6— 4,3  scheme  is  used  in  the  wall -normal  direction. 
Therefore,  the  boundary  values  are  formally  fifth-order  accurate  in  the  streamwise  direction 
since  the  points  are  equally  spaced.  However,  only  third-order  accuracy  is  maintained  in  the 
wall-normal  direction  where  the  grid  points  are  tightly  clustered  or  highly  stretched. 

The  spatial  evolution  of  the  transitional  and  turbulent  fiowfield  is  simulated  since  this  is 
the  only  physically  representative  condition  for  complex  geometries.  Inflow-outflow  boundary 
conditions  are  employed.  At  the  inflow  boundary,  the  flow  is  assumed  to  be  a  uniform 
freestream  flow  is  the  leading  edge  of  the  body  is  simulated;  otherwise  a  boundary  layer 
profile  is  used  at  the  inflow  boundary.  At  the  outflow  and  outer  boundaries,  characteristic 
boundary  conditions  axe  employed. 
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In  the  spanwise  direction,  periodic  boundary  conditions  are  employed  in  the  finite  dif¬ 
ference  or  compact  difference  framework.  At  the  wall,  no-slip  conditions  are  imposed,  the 
temperature  is  calculated  from  an  adiabatic  wall  condition,  and  either  a  zero-normal  gradient 
in  the  pressure  is  assumed  or  continuity  is  imposed  at  the  wall.  The  equation  of  state  is  also 
employed  to  calculate  either  the  density  or  pressure. 

A  fourth-order  artificial  dissipation  operator  has  been  added  to  the  second-order  finite- 
volume  formulation  to  damp  spurious  oscillations.  A  spatial  filter  was  developed  for  near¬ 
boundary  damping.  No  dissipation  operator  has  yet  been  added  to  the  higher-order  compact 
differencing  scheme,  but  a  sixth-order  filtering  is  planned  so  that  reasonable  grids  can  be 
used. 

Disturbances  can  be  introduced  into  the  computational  domain  either  through  a  suc¬ 
tion/blowing  strip  on  the  wall  or  as  a  distributed  function  at  the  inflow  boundary. 

4  NUMERICAL  VALIDATIONS 
4.1  Burger’s  Equation 

Both  the  fourth-order  compact  differences  and  second-order  finite  differences  for  the  spa¬ 
tial  derivatives  have  been  tested  for  a  model  problem.  A  comparison  of  fourth-order  com¬ 
pact  differencing  techniques  with  second-order  finite  difference  techniques  was  performed  for 
Burger’s  equation.  It  was  found  that  the  boundary  condition  treatment  has  a  significant  im¬ 
pact  on  the  order  of  accuracy  of  the  numerical  solution.  The  behavior  of  compact  difference 
techniques  is  similar  to  that  of  spectral  methods  in  that  a  larger  region  of  the  flowfield  is 
affected  by  the  boundary  condition  treatment  and  the  error  distribution  is  more  global.  The 
errors  are  more  local  to  regions  of  high  gradients  in  the  finite- difference  solution.  In  addition, 
round-off  errors  can  become  larger  than  the  truncation  error  for  small  grids  and  thus  reduce 
the  effectiveness  of  the  higher-order  differencing.  As  compact  differencing  techniques  have 
additional  CPU  requirements,  there  may  not  be  a  payoff  for  a  very  fine  grid.  The  results 
from  the  analysis  of  Burger’s  equation  have  provided  useful  information  for  the  complete 
Navier-Stokes  equations  and  highlights  of  some  of  the  test  cases  are  given  below. 

The  first  analysis  is  for  the  linear,  steady,  viscous  solution  of  Burger’s  equation  with 
Dirichlet  boundary  conditions  imposed: 

du  1  d2u 

C~di  =  lle'di 2  (21) 

«(0)  =  1,  «(1)  =  0,  c  =  1,  Re  =  10. 

Figure  1  shows  the  rms  error  in  the  calculated  solution  when  compared  with  the  exact 
solution  of  Burger’s  equation  using  standard  second-order  central  differences  and  fourth- 
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order  compact  differences.  The  slope  of  the  curves  verify  the  solution  accuracy,  except  at 
the  smaller  step  sizes  when  the  round-off  error  overtakes  the  truncation  error.  Figure  2 
shows  the  spatial  distribution  of  the  error  for  four  different  size  uniform  grids.  The  error 
in  the  second-order  finite  difference  solution  is  very  large  at  the  boundary,  but  decreases 
for  the  smaller  step  sizes  as  expected  (imax  indicates  the  number  of  grid  points  used  in 
the  computation).  The  error  in  the  fourth-order  compact  difference  solution  (with  a  fourth- 
order  one-sided  difference  for  the  first  derivative  and  a  third-order  one-sided  difference  for 
the  second  derivative  at  the  boundary)  does  not  change  as  the  grid  size  is  reduced  since  the 
round-off  error  is  higher  than  the  truncation  error  for  these  smaller  grids. 

A  second  analysis  was  performed  to  determine  the  effects  of  periodic  boundary  conditions 
for  the  unsteady,  viscous  Burger’s  equation  (c  =  1  and  Re  =  6283).  The  error  analysis  is 
shown  in  Figure  3  for  both  the  finite  difference  and  compact  difference  solutions.  The  results 
display  that  second-  and  fourth-order  accuracy  is  achieved.  Unlike  the  previous  case,  no  error 
is  introduced  at  the  boundaries.  As  periodic  boundary  conditions  are  being  employed  in  the 
spanwise  direction  in  the  full  simulation  code,  this  study  lends  confidence  to  the  capabilities 
of  the  differencing  techniques  employed. 

Another  analysis  of  the  steady,  nonlinear,  viscous  Burger’s  equation,  which  more  closely 
models  the  characteristics  of  the  Navier-Stokes  equations,  was  performed.  The  same  con¬ 
ditions  were  used  as  presented  in  the  first  example,  but  with  c  =  u.  The  error  is  shown 
in  Figure  4.  The  finite  difference  solution  shows  second-order  accuracy,  while  the  compact 
difference  solution  only  maintains  a  3.5  order  accurate  solution  as  the  step  size  is  reduced. 
The  accuracy  of  the  compact  difference  solution  was  degraded  from  fourth-order  because  the 
the  third-order  boundary  condition  for  the  second  derivative  influenced  the  flow  more  for  the 
nonlinear  Burger’s  equation.  For  very  coarse  grids,  the  accuracy  of  the  compact  difference 
solution  is  no  better  than  the  second-order  finite  difference  solution. 

The  final  example  of  the  analysis  of  Burger’s  equation  is  for  the  unsteady,  combined 
nonlinear,  viscous  case: 


«(M)  =  1, 


du 

~di 


+  (c  +  bu) 


du 

dx 


1  d2u 
Re  dx 2 


u(l,t)  =  0,  c  =  0.5,  b  =  — 1, 


Re  =  10 . 


(22) 


The  solution  of  this  equation  is  shown  in  Figure  5  and  the  error  distribution  is  shown  in 
Figure  6.  Of  particular  interest  is  the  more  localized  concentration  of  error  for  the  finite 
difference  solution,  whereas  the  errors  in  the  compact  difference  solution  are  spread  over  a 
larger  portion  of  the  domain,  similar  to  the  error  distributions  seen  in  spectral  methods.  A 
grid  density  study  showed  that  second-order  accuracy  was  achieved  for  the  finite  difference 
solution  and  fourth-order  accuracy  for  the  compact  difference  solution.  This  test  case  lends 
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confidence  to  the  higher-order  differencing  techniques  used  in  the  calculation  of  flows  with 
large  gradients,  as  will  be  present  in  the  full  simulations  of  supersonic  turbulent  flow  over  a 
curved  surface. 

The  results  presented  for  Burger’s  equation  were  obtained  on  a  uniform  grid.  However, 
the  grid  can  play  a  large  role  in  the  overall  accuracy  of  the  simulation.  Special  care  is  taken 
to  preserve  the  higher-order  accuracy  of  the  numerical  scheme.  The  grid  generation  for  the 
generalized  code  is  accomplished  through  the  MACGS  grid  generation  system  developed 
by  the  Computational  Fluid  Dynamics  Group  at  McDonnell  Douglas  Aerospace.  The  grid 
stretching  is  performed  at  a  small,  constant  rate,  as  this  yields  smaller  errors  than  if  a  tanh 
stretching  were  used.  As  Cain  (1992)  has  pointed  out,  when  using  second-  and  higher-order 
schemes  to  calculate  the  metrics,  the  metrics  will  be  exact  on  this  type  of  grid.  Also,  Cain 
has  determined  maximum  stretching  rates  based  on  the  one- dimensional  wave  equation.  If 
the  stretching  is  too  large,  the  waves  can  propagate  in  the  wrong  direction.  The  metric 
calculation  is  implemented  to  be  consistent  with  the  order  of  accuracy  of  the  flow  solver 
scheme. 

4.2  Mean  Flow 

Direct  or  large-eddy  simulation  methods  require  very  accurate  resolution  of  the  underlying 
mean  flow  to  capture  the  physics  of  transitional  or  turbulent  flows.  Great  care  has  been  taken 
to  insure  that  the  laminar  mean  flow  solutions  are  of  high  accuracy.  Comparison  has  been 
made  between  the  results  of  the  code  developed  in  this  work  and  the  compressible  similarity 
solution.  The  computational  domain  did  not  contain  the  leading  edge  of  the  plate  and  the 
compressible  similarity  solution  was  used  as  an  inflow  condition  and  to  initialize  the  flowfield. 
A  comparison  between  the  compressible  similarity  solution  and  the  newly  developed  code 
for  a  Afoo  =  3  supersonic  boundary  layer  is  shown  in  Figure  7  for  the  streamwise  u  velocity, 
the  normal  v  velocity,  and  the  temperature  T.  Excellent  agreement  is  obtained,  particularly 
for  the  normal  v  velocity. 


5  SIMULATION  RESULTS 

5.1  Large-Eddy  Simulation  of  Supersonic  Boundary  Layers 

An  investigation  of  SGS  models  for  use  in  compressible  wall  bounded  flows  was  performed 
using  a  temporal  LES  code,  CMPTBL  (CoMPressible,  Temporal,  Boundary  Layer),  devel¬ 
oped  by  Pruett  and  Zang  (1992).  This  code  was  used  to  perform  the  first  LES  of  supersonic, 
wall-bounded  turbulent  flow  (Krai  and  Zang,  1992).  However,  there  several  outstanding 
issues  concerning  SGS  models  for  supersonic  boundary  layers  arose.  This  work  was  therefore 
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continued  in  a  cooperative  effort  with  Dr.  Thomas  A.  Zang  at  NASA  Langley  Research  Cen¬ 
ter.  As  the  temporal  code  requires  much  less  CPU  time  (due  to  the  streamwise  periodicity 
approximation)  it  provided  a  good  framework  to  examine  key  numerical  issues  while  the  fully 
spatial  code  was  under  development.  The  subgrid-scale  viscosity  can  behave  similar  to  an 
artificial  viscosity  operator  for  central  differencing  codes.  As  there  is  no  artificial  viscosity 
or  upwinding  technique  in  CMPTBL,  the  values  of  the  compressible  constant  in  the  SGS 
models  were  set  high  to  provide  damping  of  unphysical  oscillations.  Since  the  earlier  work  by 
Krai  and  Zang  (1992),  Zang  (1992)  has  made  several  modifications  that  improved  stability 
and  provided  for  more  physical  damping.  These  modifications  included  treating  the  viscous 
terms  directly  as  second  derivatives  instead  of  products  of  first  derivatives.  In  the  previous 
work,  the  differencing  of  the  product  of  the  first  derivatives  did  not  damp  the  high  frequency 
waves.  This  modification  provided  for  a  much  more  stable  algorithm  for  temporal  simula¬ 
tions  of  transition  with  no  affect  on  the  solution.  In  addition,  boundary  condition  upgrades 
for  the  sixth-order  compact  differences  used  in  the  wall  normal  direction  were  made. 

Therefore,  the  earlier  findings  on  the  SGS  models  were  re-examined,  especially  the  in¬ 
fluence  of  the  compressible  constant.  Two  SGS  models  with  compressible  formulations 
were  evaluated.  One  model  is  a  compressible  formulation  of  the  Smagorinsky  model  due 
to  Speziale  et  a 1.  (1988)  and  the  second  model  is  the  structure  function  model  of  Comte  et 
aJ.  (1990).  Significant  differences  were  observed  in  the  magnitude  and  shape  of  the  fluctuat¬ 
ing  components  using  the  two  different  SGS  models.  Figure  8  shows  a  comparison  between 
results  of  a  LES  with  the  compressible  Smagorinsky  SGS  model  at  M0 0  =  4.5  and  Reg  =  5400 
using  the  earlier  version  of  the  code  (labeled  8/91  on  the  plot)  and  the  current  version  (la¬ 
beled  7/92  on  the  plot).  Also  shown  is  a  comparison  with  the  experimental  data  of  Coles 
(1953)  at  this  Mach  number  and  Reg.  All  simulations  presented  here  were  run  with  32 
Fourier  modes  in  the  streamwise  direction,  18  Fourier  modes  in  the  half-span  (with  spanwise 
symmetry  assumed),  and  72  points  in  the  wall-normal  direction  using  sixth-order  compact 
differencing.  The  same  Smagorinsky  constant  was  used  in  both  cases,  =  0.15.  Differences 
in  the  two  simulations  axe  most  evident  in  the  near  wall  region  for  the  first  derivative  of  the 
mean  streamwise  velocity  and  for  u'rm8.  The  differences  between  these  two  cases  were  not 
expected  to  be  significant.  The  intention  was  to  attempt  to  lower  Cr  with  the  more  stable 
numerical  scheme  implemented  for  the  second  derivatives,  as  the  constant  was  set  high  due 
to  the  need  for  artificial  dissipation.  With  the  modifications,  was  lowered  by  50%  to  0.075. 
Thus  the  need  for  artificial  dissipation  is  not  as  high  in  the  current  version  of  CMPTBL. 
Results  with  the  lower  values  of  Cr  are  shown  in  Figure  9.  The  slope  of  the  mean  velocity  and 
temperature  changes  significantly  at  the  wall.  With  the  decrease  in  Cr,  the  SGS  viscosity  is 
reduced,  which  leads  to  higher  values  for  u'rmt.  Similar  results  were  found  for  the  structure 
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function  SGS  model  and  results  are  shown  in  Figure  10.  The  constant  was  lowered  by 
30%  from  0.30  to  0.21  before  numerical  stability  became  an  issue. 

As  the  proper  value  for  this  constant  is  not  known  for  supersonic,  wall-bounded  turbu¬ 
lent  flows,  the  results  of  these  temporal  simulations  can  be  used  as  a  guide  in  the  spatial 
simulations. 

Figure  11  shows  results  from  the  temporal  LES  of  supersonic  turbulent  flow  over  a  flat 
plate  at  =  4.5.  Shown  is  a  constant  density  surface  at  one  instant  in  the  simulation 
using  the  compressible  Smagorinsky  SGS  model.  The  simulation  appears  to  be  modeling 
the  important  physics  of  the  boundary  layer  as  evidenced  by  the  similarity  of  certain  flow 
features  to  those  observed  in  experiments.  Horseshoe-shaped  structures,  visible  in  the  fig¬ 
ure,  are  almost  identical  to  structures  found  in  incompressible  boundary  layers  (see,  e.g. 
Robinson,  1991)  and  are  very  similar  to  those  observed  in  compressible  flows  where  much 
less  experimental  information  is  available  (Spina  et  a/.,  1991).  Additional  visualization  of 
low  momentum  fluid  near  the  wall,  as  shown  in  Figure  12  indicates  low-speed  streaks  and 
shows  regions  of  these  streaks  lifting  away  from  the  wall  in  the  vicinity  of  the  horseshoe¬ 
like  structures,  as  observed  experimentally.  In  Figure  13,  contours  of  density  at  a  spanwise 
plane  in  the  boundary  layer  are  shown,  convecting  with  0.85Z7e  which  is  the  approximate 
convection  speed  of  large-scale  motions.  Superimposed  are  velocity  vectors.  Low  density 
fluid  is  lifted  up  from  the  wall  and  local  shear  layers  are  observed  away  from  the  wall  that 
are  associated  with  the  low-density  structures.  The  downstream  “bulge”  is  in  early  stages  of 
development  as  depicted  by  the  strong  shear  layers  seen  on  the  upstream  side.  The  results 
are  in  qualitative  agreement  with  incompressible  motions. 


5.2  Direct  Numerical  Simulation  of  Gortler  Flow 

Counter-rotating  streamwise  vortices  are  observed  in  transitional  and  turbulent  bound¬ 
ary  layers  and  play  an  important  role  in  the  dynamics  governing  the  dynamics  governing 
the  growth,  breakdown,  and  transition  to  turbulence.  The  streamwise  vortices  develop  in 
boundary-layer  flow  over  a  concave  wall  due  to  the  imbalance  of  the  centrifugal  forces  and 
the  radial  pressure  gradient  as  first  predicted  by  Gortler  (1940).  As  this  Gortler  instability 
is  a  common  phenomena  in  various  flows  involving  curved  boundaries,  it  was  chosen  to  be 
our  validation  case  for  the  numerical  method  developed  to  date. 

The  geometry  and  flow  conditions  closely  approximate  the  incompressible  experiments  of 
Swearingen  and  Blackwelder  (1987).  A  sketch  of  the  streamwise  vortices  developing  on  a 
concave  wall  is  shown  in  Figure  14.  The  spanwise  wavelength  of  the  vortices  is  denoted  by 
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A  in  the  figure.  The  governing  nondimensional  parameter  is  Gortler  number  defined  by 

where  6  is  the  momentum  thickness,  is  the  freestream  velocity,  R  is  the  radius  of  cur¬ 
vature,  and  v  is  the  kinematic  viscosity.  In  contrast  to  the  experiments,  the  simulations 
were  performed  at  =  0.5  in  order  to  validate  the  code  at  subsonic  Mach  numbers.  The 
radius  of  curvature  is  R  —  3.2  m,  the  freestream  velocity  is  Uoo  —  174  m/a,  the  freestream 
temperature  is  T  =  300  K,  the  unit  Reynolds  number  per  meter  is  Re  =  3.2  x  10s  m  *, 
and  the  streamwise  length  of  the  computational  domain  is  L  =  251  cm.  The  span  wise  wave¬ 
length  is  chosen  to  match  the  experimentally  observed  wavelength,  A  =  2.3  cm.  The  flow 
is  unstable  for  Gog  >  0.3  and  the  critical  value  for  transition  to  turbulence  occurs  between 
6  <  Gog  <  10  or  78  cm  <  x  <  152  cm.  The  disturbance  is  introduced  into  the  computational 
domain  at  the  inflow  boundary  through  the  introduction  of  perturbations  in  the  spanwise 
and  wall-normal  velocities. 

The  grid  for  this  direct  simulation  was  generated  using  MACGS,  with  constant  stretching 
in  the  streamwise  direction  and  constant  spacing  in  the  spanwise  direction.  A  small,  con¬ 
stant  rate  stretching  is  used  in  the  wall  normal  direction.  The  computational  grid  is  shown  in 
Figure  15.  The  grid  is  104  X  61  x  31.  Curvature  starts  at  the  leading  edge  and  the  computa¬ 
tional  domain  extends  from  5  cm  upstream  of  the  leading  edge  to  251  cm  downstream  of  the 
leading  edge  along  the  surface.  The  spanwise  extent  of  the  domain  is  2.3  cm  corresponding 
to  one  wavelength. 

The  initial  stages  of  this  instability  exhibit  the  expected  exponential  growth  in  the  stream- 
wise  direction.  Figure  16  shows  the  downstream  growth  of  the  amplitude  of  the  spanwise 
velocity  near  the  wall  at  a  value  of  y+  «  10.  The  disturbance  has  propagated  to  110  cm  and 
the  flowfield  downstream  of  this  location  is  not  yet  developed. 

Normal  profiles  of  the  mean  streamwise  velocity  at  five  downstream  locations  are  shown 
in  Figure  17  for  spanwise  locations  corresponding  to  a  peak  and  a  valley  region.  A  ‘peak’  is 
the  region  where  u'(z )  is  a  maximum  and  U(z)  is  a  minimum.  A  noticeable  distortion  of  the 
mean  laminar  flow  is  observed  in  the  figure.  The  peak  region  shows  a  velocity  defect  and 
the  valley  region  shows  an  increase  in  the  mean  velocity.  Concentrated  regions  of  low-  and 
high-speed  fluid  are  produced  and  the  boundary  layer  has  a  wavy  appearance  in  the  spanwise 
direction.  In  the  low-speed  region  at  x  =  80  cm,  the  mean  velocity  profile  is  inflectional  in 
the  normal  y  direction.  The  inflectional  profiles  indicate  the  imminent  onset  of  a  tertiary 
instability  and  the  breakdown  to  turbulence.  The  mean  streamwise  profiles  are  in  agreement 
with  the  experiments  of  Swearingen  and  Blackwelder  (1987),  in  which  the  inflectional  profile 
was  observed  also  at  x  =  80  cm. 
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The  low-  and  high-speed  regions  in  the  spanwise  direction  results  in  an  alternating  increase 
and  decrease  in  the  thickness  of  the  boundary  layer.  The  calculated  development  of  the 
displacement  thickness  8*  in  the  downstream  direction  for  the  low-  and  high-speed  regions 
is  compared  with  the  laminar  growth  in  Figure  18.  Also  shown  are  the  measurements  of 
Swearingen  and  Blackwelder.  In  the  low-speed  region,  8*  increases  dramatically  near  x  — 
40  cm.  A  corresponding  decrease  in  8 *  is  observed  in  the  high-speed  region.  Again,  the 
Gortler  instability  has  reached  about  x  —  110  cm  in  the  simulations.  Overall,  the  agreement 
between  the  calculations  and  the  experiments  is  quite  good. 

The  growth  of  the  vortex  system  is  illustrated  in  Figure  19.  Iso-contours  of  the  streamwise 
velocity  are  shown  at  several  downstream  locations.  Also  shown  is  a  comparison  with  the 
experiments  of  Swearingen  and  Blackwelder  at  the  same  locations.  Initially  the  disturbances 
are  linear  and  the  boundary  layer  appears  unperturbed,  but  is  quickly  modified  downstream, 
as  the  disturbance  becomes  non-linear  and  the  influence  of  the  counter-rotating  vortices 
is  evident.  As  low-momentum  fluid  is  pumped  away  from  the  wall  between  vortices,  the 
boundary  layer  and  on  either  side  of  the  domain  the  boundary  layer  is  thinned  due  to  the 
induction  of  fluid  towards  the  wall.  The  comparison  between  the  direct  numerical  simulations 
and  the  experiments  is  qualitatively  similar.  However,  it  appears  that  the  development  of 
the  Gortler  vortices  in  the  simulations  is  delayed  10  —  20  cm  relative  to  the  experiments.  This 
result  is  not  unexpected  because  factors  present  in  the  experiment,  such  as  surface  roughness 
and  possibly  larger  amplitude  disturbances,  will  cause  the  boundary  layer  to  be  less  stable. 

Streamwise  vortices  initiated  by  the  Gortler  instability  exist  in  the  presence  of  a  strong 
mean  gradient  of  the  streamwise  velocity  U(y).  Between  the  two  vortices,  the  induced 
motion  removes  low-momentum  fluid  away  from  the  wall.  At  one-half  wavelength  away  in 
the  spanwise  direction,  the  vortices  pump  high-speed  fluid  towards  the  wall.  The  counter¬ 
rotating  vortices  are  illustrated  in  Figure  20,  which  shows  contours  of  spanwise  velocity.  The 
low-speed  regions  that  result  from  this  lift-up  of  low-momentum  fluid  develop  inflectional 
velocity  profiles  as  shown  in  Figure  17. 

Although  these  simulations  have  not  yet  reached  a  fully  turbulent  state,  at  the  time  of 
this  report  the  Gortler  instability  is  still  proceeding  downstream.  It  is  anticipated  that  the 
final  breakdown  to  turbulence  will  occur. 

6  SUMMARY 

The  objective  of  this  work  has  been  to  develop  a  DNS  capability  for  compressible,  non- 
canonical  wall-bounded  flows  which  is  directly  applicable  to  actual  flight  vehicles.  A  signif¬ 
icant  part  of  this  effort  was  placed  on  the  development  of  the  numerical  algorithm  that  is 
used  to  perform  DNS  of  compressible,  turbulent  flows  over  bodies  with  surface  curvature. 
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Much  of  the  work  has  therefore  focused  on  resolving  issues  necessary  in  formulating  the 
numerical  method.  A  second-order  finite  volume  approach  and  both  fourth-order  and  sixth- 
order  compact  differences  are  available  for  the  spatial  derivatives.  A  third-order  low-storage 
Runge-Kutta  scheme  is  used  for  the  time  integration.  The  algorithm  allows  for  general¬ 
ized  coordinates  so  that  simulations  about  complex  geometries  can  be  performed.  Extensive 
testing  was  performed  for  two  SGS  models,  the  compressible  Smagorinsky  and  structure 
function  models,  to  determine  the  effects  of  the  SGS  constant.  The  study  reinforced  the 
need  for  high-order  dissipation  models  in  the  simulation  code,  as  the  SGS  viscosity  also 
behaves  as  an  artificial  viscosity  for  central  differencing  codes.  The  numerical  method  was 
validated  for  flow  over  a  concave  surface  at  =  0.5  in  which  streamwise  vortices  develop 
in  the  boundary-layer  due  to  the  Gortler  instability.  The  geometry  and  flow  conditions 
closely  approximated  the  experiments  of  Swearingen  and  Blackwelder  (1987).  The  simula¬ 
tions  captured  the  essential  features  of  the  experiments  in  which  counter-rotating  vortices 
developed  near  the  wall.  The  laminar  boundary  layer  on  the  concave  wall  rapidly  becomes 
three-dimensional  as  these  streamwise  vortices  develop.  The  vortices  grow  downstream  and 
higher-momentum  fluid  is  pulled  towards  the  wall  from  the  outer  flow.  The  low-speed  fluid 
lying  between  the  vortices  induces  low-momentum  fluid  away  from  the  wall,  leading  to  in¬ 
flectional  streamwise  velocity  profiles,  in  good  agreement  with  the  incompressible  studies  of 
Swearingen  and  Blackwelder. 
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Figure  1:  Rms  error  in  computing  the  steady,  linear,  viscous  Burger’s  equation  at  Re  ■■ 
second-order  finite  differences  and  fourth-order  compact  differences. 
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Figure  2:  Spatial  distribution  of  error  for  the  steady,  linear,  viscous  Burger’s  equation  at 
using  second-order  finite  differences  and  fourth-order  compact  differences  for  four  different 
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Figure  3:  Rms  error  in  computing  the  unsteady,  linear,  viscous  Burger’s  equation  at  Re  =  6283  and  c  =  1 
with  periodic  boundary  conditions  using  second-order  finite  differences  and  fourth-order  compact  differences. 


Figure  4:  Error  in  computmg  the  unsteady,  nonlinear,  viscous  Burger’s  equation  at  Re  =  10  and  tt  =  c  using 
second-order  finite  differences  and  fourth-order  compact  differences. 
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Figure  6:  Spatial  distribution  of  error  for  the  unsteady,  nonlinear,  viscous  Burger’s  equation  at  Re  —  10, 
c  =  0.5,  and  b  =  -1  using  second-order  finite  differences  and  fourth-order  compact  differences  and  Ax  =  O.Ol! 
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Figure  7:  Comparison  of  laminar  mean  flow  solution  with  the  compressible  similarity  solution  for  a  =  3 
flat  plate  boundary  layer. 


Figure  8:  Influence  of  modifications  in  the  high-order  differencing  scheme  on  the  mean  flow,  the  SGS  viscosity, 
and  the  rms  of  the  streamwise  velocity  using  the  compressible  Smagorinsky  SGS  model  (cr  =  0.15)  for  a 
temporal  LES  of  a  M «,  =  4.5  turbulent  boundary  layer  at  Reg  =  5400.  Comparison  is  shown  with  the 
experimental  data  of  Coles  (1953). 
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Figure  9:  Influence  of  the  compressible  Smagorinsky  constant  on  the  mean  and  fluctuating  flow  of  a  supersonic 
turbulent  boundary  layer  using  the  modified  temporal  LES  code  with  the  compressible  Smagorinsky  SGS 
model.  Comparison  is  also  shown  with  the  experimental  data  of  Coles  (1953). 
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Figure  10:  Influence  of  the  compressible  SGS  constant  on  the  mean  and  fluctuating  flow  of  a  supersonic 
turbulent  boundary  layer  using  the  modified  temporal  LES  code  with  the  compressible  structure  function 
SGS  model.  Comparison  is  also  shown  with  the  experimental  data  of  Coles  (1953). 


Figure  12:  Low  momentum  (streamwise)  fluid  near  the  wall  from  a  temporal  LBS  of  a  supersonic  turbulent 
boundary  layer  using  the  compressible  Smagorinsky  SGS  model.  Low-speed  streaks  are  seen  lifting  away 
from  the  wall. 
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Figure  13:  Velocity  vectors  and  density  contours  identifying  structural  features  in  a  Mach  4.5  turbulent 
boundary  layer. 


Figure  14:  Sketch  of  the  streamwise  vortices  developing  on  a  concave  wall  due  to  the  Gortler  instability 
mechanism  (from  Swearingen  and  Blackwelder  (1987)). 
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Figure  18:  Downstream  development  of  the  displacement  thic 
Also  shown  are  the  experimental  measurements  of  Swearingei 


30 


MDC  94X0041 


x  =  90  cm 


Z  Z  m  I’-f- 


x  =  1 00  cm 


x  =  1 1 0  cm 


Figure  19:  Iso-contours  of  the  mean  streamwise  velocity  in  the  cross-stream  (y,  z)-plane 
for  several  streamwise  x  locations.  Also  shown  on  the  right  side  are  the  experimentally 
determined  contours  of  Swearingen  and  Blackwelder  (1987). 


